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METHOD AND SYSTEM FOR MARKER LOCALIZATION 



CROSS-REFERENCE TO RELATED APPLICATIONS 

[0001] This application is related to U.S. Patent Application No. 10/334,700 filed 

December 30, 2002, Attorney Docket No. 341 14-8003.US00, and U.S. Patent 
Application No. 10/382,123, filed March 4, 2003, Attorney Docket No. 34114- 
8008. US01, both of which are incorporated herein by reference in their entirety. 

BACKGROUND 

[0002] Implantable markers have been used to identify locations within objects, 

such as a human body. For example, a marker may be implanted in a patient 
within an organ of interest. As the patient moves, the marker can be used to track 
the location of the organ. Various techniques have been used to identify the 
location of such markers. For example, one technique requires a person to move 
a sensor over the area above the marker. When the sensor is positioned directly 
over the marker, the person may be given a visual or audio indication. A difficulty 
with such a technique for identifying a marker location is it does not provide the 
actual location of the marker (e.g., x, y, and z coordinates) and may not have the 
needed degree of accuracy. 

[0003] One technique has been described for locating a source of an electrical 

signal within the brain by sensing the magnetic field generated by an electrical 
signal. This technique relies on a probability distribution to help identify the 
location of the source. A difficulty with such a technique is it may not have the 
needed degree of accuracy. 

[0004] It would be desirable to have a technique for locating an implantable marker 

in real time and in the needed degree of accuracy. 
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BRIEF DESCRIPTION OF THE DRAWINGS 



[0005] Figure 1 is a block diagram illustrating a bounding volume containing a 

marker in one embodiment. 
[0006] Figure 2 is a block diagram of a bounding box that is divided into a cubic 

grid of four sub-boxes in each of the x, y, and z directions. 
[0007] Figure 3 illustrates components of the location system in one embodiment. 

[0008] Figure 4 is a flow diagram illustrating the processing of the find bounding 

sub-box component in one embodiment. 
[0009] Figure 5 is a flow diagram illustrating the processing of the find closest grid 

point component in one embodiment. 
[0010] Figure 6 is a flow diagram illustrating the processing of the find adjacent 

sub-box component in one embodiment. 
[0011] Figure 7 is a flow diagram illustrating the processing of the find neighbor 

maximum point component in one embodiment. 
[0012] Figure 8 illustrates an interpolation point within a bounding sub-box. 

[0013] Figure 9 is a flow diagram illustrating the processing of the find marker point 

within a bounding sub-box component in one embodiment. 

DETAILED DESCRIPTION 

[0014] A method and system for locating a source of a signal is provided. In one 

embodiment, the source is a wireless marker that is implanted in an object, such 
as a human body. The signal generated by the marker is a magnetic field. The 
location system uses an array of sensors to measure the magnetic field of the 
marker at various sensor locations. Each sensor generates a measurement of 
one component of the magnetic field integrated over the aperture of the sensor. 
The location system determines the location of the marker (i.e., marker location) 
from a set or array of measurements taken from the sensors (i.e., set of actual 
measurements). The location system compares the set of actual measurements 
to sets of reference measurements for various known locations within a bounding 
volume (also referred to as a localization volume). The bounding volume delimits 
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the three-dimensional area in which the marker can be localized. A reference 
measurement for a known location indicates the measurements to be expected 
from the sensors when the marker is located at that known location. Based on the 
comparisons, the location system identifies the set of reference measurements 
that most closely matches the set of actual measurements. The known location of 
the identified set of reference measurements represents the known location that is 
closest to the marker location, which is referred to as the "closest known location." 
The location system then uses sets of reference measurements for known 
locations near the closest known location to more accurately determine the marker 
location when it is not actually at one of the known locations. In one embodiment, 
the location system determines the marker location based on an interpolation of a 
set of calculated measurements from the sets of reference measurements of 
known locations near the closest known location. Thus, the location system uses 
the set of reference measurements to find a known location that is close to the 
marker location to an accuracy that is dependent on the spacing of the known 
locations. The location system then uses an interpolation of sets of reference 
measurements at known locations near the closest known location to more 
accurately identify the marker location at a location between the known locations. 
[0015] Figure 1 is a block diagram illustrating a bounding volume containing a 

marker in one embodiment. The bounding volume in this embodiment is a 
rectangular volume referred to as a bounding box. Locations or points within the 
bounding box 100 are specified using x, y, and z coordinates relative to the origin 
101 with coordinates (0,0,0). The bounding box has X units along the x axis, Y 
units along the y axis, and Z units along the z axis. The marker location 102 has 
coordinates (Ax, Ay, Az). The location system divides the bounding box into a grid 
of sub-boxes of a uniform size. Figure 2 is a block diagram of a bounding box that 
is divided into a cubic grid of four sub-boxes in each of the x, y, and z directions. 
The bounding box 200 thus contains 64 sub-boxes. Although the bounding box 
and the sub-boxes shown in this embodiment are cubes, one skilled in the art will 
appreciate that bounding boxes and sub-boxes of various shapes can be used. 
Sub-box 201 has grid points (0,0,0), (0,0,1), (0,1,0), (0,1,1), (1,0,0), (1,0,1), (1,1,0), 
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and (1 ,1 ,1 ) at its corners. Each corner point of a sub-box corresponds to a known 
location and is referred to as a "grid point." Each grid point, or more generally a 
reference point, thus has a set of reference measurements associated with it. The 
sub-boxes are identified by the corner point with the lowest x, y, and z coordinates 
that conceptually corresponds to the lower left forward corner of a sub-box. Sub- 
box 201 is uniquely identified by grid point (0,0,0), and the sub-box next to it to the 
right is uniquely identified by grid point (1,0,0). 
[0016] The location system uses an objective function to evaluate which of the sets 

of reference measurements is most consistent with a set of actual measurements. 
Using the objective function, the system can locate the closest grid point to the 
marker point (i.e., the point corresponding to the marker location) using various 
search techniques. One search technique is a comprehensive search that 
calculates the objective function at each grid point and selects the grid point 
whose set of reference measurements is most consistent with the set of actual 
measurements based on an evaluation of the objective function. Another search 
technique starts at a grid point and evaluates the objective function at the grid 
point and at neighbor grid points. The search technique then selects the neighbor 
grid point whose objective function evaluation indicates it is closer to the marker 
location and repeats the process of evaluating the objective function of its 
neighbor grid points. This process is repeated until a grid point is found whose 
objective function evaluation indicates that its set of reference measurements is 
more consistent than any of its neighbor grid points to the actual measurements 
(e.g., closer to the marker point). Another search technique divides the bounding 
volume into sub-volumes and evaluates the objective function at the center grid 
point of each sub-volume. The search technique selects the sub-volume whose 
objective function evaluation for its center grid point indicates that it is closer to the 
marker point than the other center grid points. The search technique then divides 
the selected sub-volume into sub-sub-volumes and repeats the evaluation of the 
objective function and selection of sub-volumes until a sub-volume cannot be 
further divided. One of the grid points of the final sub-volume is the closest grid 
point. 
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[0017] In one embodiment, once the closest grid point is located, the location 

system identifies which of the eight sub-boxes that include the closest grid point as 
a corner contains the marker location. This sub-box is referred to as the 
"bounding sub-box." (More generally, the location system identifies a sub-volume 
near the closest reference point that presumably contains the marker point as the 
bounding sub-volume.) In one embodiment, the location system evaluates the 
objective function at the two adjacent grid points to the closest grid point in each of 
the x, y, and z directions. The location system selects the adjacent grid point in 
each direction whose evaluation of the objective function indicates it is closer to 
the marker point. The location system then selects the sub-box that has all three 
of the selected adjacent grid points as corners of the sub-box that contains the 
marker point. For example, if coordinates of the closest grid point are (10, 10, 10), 
the adjacent grid points in the x direction are (9, 10, 10) and (11, 10, 10), in the y 
direction are (10, 9, 10) and (10, 11, 10), and in the z direction are (10, 10, 9) and 
(10, 10, 11). If the closer adjacent grid points as indicated by the objective 
function are (9, 10, 10) for the x direction, (10, 11, 10) for the y direction, and (10, 
10, 11) for the z direction, then the location system selects the sub-box whose 
lower left forward corner is at grid point (9, 10, 10) as the bounding sub-box. 

[0018] Once the location system has identified the bounding sub-box, it then 

searches for the marker point within the bounding sub-box. Since the location 
system does not have a set of reference measurements for the points within the 
bounding sub-box, except for the corner points, the location system cannot 
evaluate the objective function directly from the sets of reference measurements 
for those points. In one embodiment, the location system can interpolate the value 
of an objective function for points within the bounding sub-box from the values of 
the objective functions at the corner points (and possibly additional grid points 
near the bounding sub-box). For example, the value of the objective function for a 
point within the bounding sub-box can be an interpolation based on the distance 
from the point to each of the corners (e.g., a trilinear interpolation of the values of 
the objective function at the corners). It has been found, however, that 
interpolation of the objective function can be computationally expensive. In 
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another embodiment, the location system can interpolate a set of reference 
measurements for a point within the bounding sub-box based on the sets of 
reference measurements at the corner points (and possibly additional grid points 
near the bounding sub-box). Such interpolated reference measurements are 
referred to as "calculated measurements." The location system uses the set of 
calculated measurements for a point to evaluate the objective function for that 
point. In one embodiment, the location system uses a trilinear interpolation of the 
sets of reference measurements at the corner points adjusted by a second order 
difference of the sets of reference measurements of the grid points adjacent to the 
corner points. One skilled in the art will appreciate that many different 
interpolation techniques can be used by the location system, for example, as 
discussed in W.H. Press, et al. Numerical Recipes in C/C++: The Art of Scientific 
Computing , (Cambridge Univ. Press 2d ed. Feb. 2002), which is hereby 
incorporated by reference. To search for the marker point in one embodiment, the 
location system uses a Taylor series expansion of the objective function using 
calculated measurements to guide the search. Conceptually, the location system 
evaluates the objective function at a candidate point and uses the gradient (i.e., 
the first derivative) and the Hessian (i.e., the second derivative) of the objective 
function to identify, as the next candidate point, the point at which the gradient is 
zero. The location system calculates the value of the objective function at that 
next candidate point and repeats the process of selecting next candidate points 
until a desired value of the objective function is achieved or until the distance 
between the successive candidate points is within a desired tolerance. 
Mathematical Description of the Location System 
[0019] The location system is provided with a vector s b containing a set of actual 

measurements arising from a marker, emitting a dipolar magnetic field, located at 
an unknown location r b and having a magnetic moment m b . The elements of 
vector s b have a one-to-one correspondence to the sensors in the array. The 
location system then estimates location r b from the actual measurements. In the 
absence of measurement error, additive noise, or other corrupting influences, the 
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relationship between the actual measurements, the marker location, and the 
magnetic moment vector of the marker is represented by 

s„=H 6 m 6 (1) 
where H 6 is a matrix containing coefficients indicative of the relative contribution of 
the x, y, and z components of the magnetic moment to each actual measurement 
when the marker is at location r b . Specifically, H t in equation (1) has three 
columns; the first column indicates the sensitivity of the measurement to the x 
component of m 6 , the second column to the y component, and the third column to 
the z component. Hence, if the location system has 32 sensors, then H fc is a 32- 
by-3 matrix. Furthermore, H 6 is sensitive to the location of the marker, but 
independent of the magnetic moment of the marker, and in general each possible 
marker location r can be associated (given a coordinate system and units of 
measurement implicit in equation (1)) with a specific matrix H. 
[0020] The H matrices of interest (corresponding to marker locations of interest) 

can be empirically or computationally determined in advance of locating a marker. 
The coefficients of H can be empirically determined for a given marker location by 
first locating a test marker at the given location and orienting it in the x direction 
and using the resulting measurements as the first column of H. The marker can 
then be oriented in the y direction and the resulting measurements used as the 
second column of H, and similarly for the z direction. It will be appreciated that 
the measurements thus obtained are scaled by the magnitude of the magnetic 
moment of the test marker; hence, this magnitude is an implicit unit of 
measurement. Alternatively, when the sensor array is comprised of coils, the 
coefficients of H can be calculated using Faraday's equation (as in, for example, 
S. Ramo, et al. Fields and Waves in Communication Electronics (John Wiley and 
Sons 2d ed. 1984), which is hereby incorporated by reference. The determination 
of H in advance of locating a marker has significant practical advantages, as an 
empirical determination may be impossible at the time of locating a marker, and a 
computational determination may be difficult or expensive to do with sufficient 
speed at the time of locating. 
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[0021] It will be appreciated that the numerical content of H determined by the 

methods above will depend on the choice of coordinate system, as well as the 
magnitude of the magnetic moment of the test marker. Specifically, a rotation of 
the Cartesian coordinates will result in an alternative H, in which the columns are 
linear combinations of the columns of H associated with nonrotated coordinates. 
These differences are unimportant to the problem of determining the location of a 
marker, arising as they do from the arbitrary imposition of a coordinate system. 
Accordingly, all H matrices associated with the same location are considered 
equivalent, even if their numerical content differs. It follows that the columns of all 
equivalent H matrices span the same subspace. See B. Noble & J.W. Daniel, 
Applied Linear Algebra (Prentice-Hall 3d ed. 1988), which is hereby incorporated 
by reference. The columns of each H matrix thus determined represent a set of 
reference measurements associated with a possible marker location, in which the 
pertinent characteristic of the set of reference measurements is the subspace thus 
spanned. In other words, in the absence of noise, interference, or other sources of 
error, an actual measurement from a marker will be a linear combination of the 
columns of the H matrix associated with the location of the marker. 

[0022] The problem of inverting equation (1), i.e., estimating the location of a 

marker given the vector of actual measurements, is the problem of finding the H 
matrix that is most consistent with the actual measurements (and any other 
relevant prior information, such as the noise or interference to which the 
measurements are subject), in which case the best estimate of marker location is 
the location associated with the most consistent H. 

[0023] The location system uses an error metric to quantify the consistency 

between a candidate H (corresponding to a candidate location r of the marker 
whose field is measured) and a vector s of actual measurements. The error metric 
is represented by the scalar 

£,=(Hm-s) r (Hm-s) (2) 
The minimization of E x with respect to unknown parameters represents a least- 
mean-square ("LMS") estimate of those parameters. See, for example, S.M. Kay, 
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Fundamentals of Statistical Signal Processing: Estimation Theory (Prentice Hall 
PTR 1993), which is hereby incorporated by reference. 
[0024] When the marker is a wireless implant, it is usually impractical to know the 

magnetic moment with sufficient precision for localization. In this case, the LMS 
solution for m in terms of H and s is well known to be 

m £M =(H r H)" 1 H r i (3) 

Substituting this into equation (2), it is seen that the minimization of E x 
corresponds to the maximization of 

£ 2 =s r H(H r H)~ ! H r s (4) 

This is the objective function used in a preferred embodiment of the invention. It 

will be appreciated that H (h 7 H B 7 is a projection matrix (see G. Strang, Linear 

Algebra and its Applications (Academic Press 2d ed. 1980)), which is hereby 

incorporated by reference), and H (h 7 H ) _1 H r s is the projection of s into the 

subspace spanned by the columns of H. Accordingly, E 2 is the inner product of s 
with its projection into the column space of H. Each equivalent H yields the same 
value of E 2 . The LMS estimate of the marker location is the location r 
corresponding to the set of reference measurements in the matrix H that effects 
the maximization of E 2 . 
[0025] In a preferred embodiment of the invention, the H matrices corresponding 

to a grid of points spanning a bounding box (i.e., a localization volume) are 
predetermined and stored in computer-readable form. At the time of locating a 
marker, a two-step process is followed. First, a coarse search is effected, in which 
equation (4) is evaluated on appropriate grid points to identify the grid point which 
maximizes E 2 , followed by the identification of eight grid points that are expected 
to bound an off-grid location which maximizes E 2 . These grid points define a 
bounding sub-box. Second, a fine search is effected, in which equation (1) is 
maximized using interpolated values of H, corresponding to off-grid locations, that 
are computed using the predetermined H matrices. 
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[0026] The inventors have found that the columns of H matrices thus interpolated 

can span the correct subspace with a high degree of accuracy, even with relatively 
coarse spatial sample grids. It has been found, conversely, that interpolated 

projection matrices H (H r H ) 1 H r do not represent subspaces with sufficient 

accuracy unless the spatial sample grid is substantially finer. The use of relatively 
coarse spatial grids has substantial practical advantages in that the amount of 
computer memory required to store the predetermined matrices is proportional to 
the number of grid points. For example, if the localization volume is 14cm x 14cm 
x 20cm, and one 32 x 3 H matrix is predetermined for each spatial grid point on a 
1cm pitch, the total storage requirement is 376,320 numbers in computer memory. 

By contrast, if one 32 x 32 H (H r H H r matrix is predetermined on a spatial grid 

with a 1mm pitch, the storage requirement exceeds four billion numbers. 

Interpolation of Reference Measurements 
[0027] In the neighborhood of any bounding sub-box, off-grid H matrices can be 

computed using any of a number of interpolation techniques known to those skilled 

in the art of numerical computation. In a preferred embodiment of the invention, a 

second order polynomial expansion is used. 
[0028] An interpolation point is denoted as a fraction (Ax,Ay,Az) away from the 

origin of the bounding sub-box. The fractions vary from (0,0,0) at the origin to 

(1,1,1) at the corner diagonal from the origin. Besides the eight points at the 

corners of the bounding sub-box, the preferred embodiment uses 24 additional 
points, which are the two points adjacent to the bounding sub-box along each of 
the 12 edges. 

[0029] Figure 8 illustrates an interpolation point within a bounding sub-box. The 

bounding sub-box has eight corners labeled H^, H^, H 010 , H ou , H 10o , H 101 , H 110 , 
H U1 . The corner labeled is at the origin. The coordinates of the interpolation 
point 801 are (Ax,Ay,Az). A trilinear interpolation combines the H matrices at the 

corners weighting them linearly based on the distances from the interpolation 
point. Thus, for example, an interpolation point in the center of the bounding sub- 
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box (e.g., (.5, 5, 5) ) can be associated with a trilinearly interpolated H matrix that 

is the average of the H matrices of the corners. More generally, a trilinear 
interpolation takes the form given by the following equation: 

H(Ax,Ay,Az) = (1 - Ax)(l - Ay)(l - Ar)H M + 

(Ax)(l-Ay)(l-Az)H 100 + 

(l-Ax)(Ay)(l-Az)H 010 + 

(l-Ax)(l-Ay)(Az)H 001 + 

(Ax)(Ay)(l-Az)H U0 + 

(Ax)(l-A^)(Az)H 

101 1 

(l-Ax)(A^)(Az)H 0I1 + 
(Ax)(A^)(Az)H in 

To add the second order terms, the location system fits a second order 
function along each of the 12 edges of the bounding sub-box using the four grid 
points along the line of each edge. For example, the second order difference 
interpolation for the edge between the corners labeled and H 100 would use 
the corners labeled H_ 100 , H,^, H 100 , and H 200 . An interpolated H matrix includes 
the eight terms of the trilinear interpolation and the 12 terms of the second order 
interpolation and is represented by the following equation: 
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H(Ax,Ay,Az) = (l - Ax)(l - Ay)(l - Az)H 000 + 
(Ax)(l-Aj,)(l-Az)H 100 + 
(l-Ax)(Ay)(l-Az)H 010 + 
(l-Ax)(l-Ay)(Az)H 001 + 
(Ax)(A^)(l-Az)H 110 + 
(Ax)(l-Ay)(Az)H 101 + 
(l-Ax)(Ay)(Az)H 0U + 
(Ax)(A^)(Az)H m + 

Ax(l - Ax)(l - Ay)(l - Az)[-H_ 100 + + H 100 - H 200 ]/4 + 
Ax(l-Ax)(A>.)(l-Az)[-H 110 +H 010 + H 110 -H 210 ]/4 + 
Ax(l-Ax)(l- A^)(Az)[-H 101 + H W1 +H 101 -H 201 ]/4 + 
Ax(l-Ax)(A.y)(Az)[-H H1 +H 0I1 +H 1U -H 211 ]/4 + 
Ay(l - Ay)(l - Ax)(l - Az)[-H^ 10 + + H 010 - H 020 ]/4 + 
AXl-Ay)(Ax)(l-Az)[-H W0 +H 100 +H 110 -H 120 ]/4 + 
Ay(\-Ay)(l- Ax)(Az)[-H 0 _ u +H 001 +H on -H 021 ]/4 + 
A^(l-Ay)(Ax)(Az)[-H 1 _ n + H 101 +H m -H 121 ]/4 + 
Az(l - Az)(l - Ax)(l - A^-H^ + H (m + H W1 - H^/4 + 
Az(l-Az)(Ax)(l-A^)[-H 1(M + H 100 + 
Az(l - Az)(l - Ax)(A^) [-H ow + H 010 + H 01 , - H 012 ]/4 + 
Az(l-Az)(Ax)(Ay)[-H lw + H 110 +H U1 -H 112 ]/4 

[0031] A characteristic of the interpolated matrices in this embodiment is that the 

values of the interpolated matrix match the set of reference measurements at the 
eight corner points of the bounding sub-box exactly, and at each of the 24 
additional points in a least-mean-squared sense. Another characteristic of this 
and similar embodiments is that the interpolated matrices can be computed very 
quickly, to effect real-time localization, on general-purpose computers without 
special hardware. A third characteristic is the simplicity of spatial derivatives of the 
interpolated matrices, which enables rapid iteration of the fine search for the 
maximum of the objective function. 

[0032] Using an interpolated matrix H associated with a point in space, the 

objective function associated with that point can be represented as: 
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[0033] 



[0034] 



E 2 = s r H (H r H )"' H r s (7) 

In a preferred embodiment, the system represents the objective function as 
a Taylor series expansion and solves for the point at which the gradient is zero. 
See R. Fletcher, Practical Methods of Optimization (John Wiley and Sons 1980), 
which is hereby incorporated by reference. A Taylor series expansion of the 
objective function about some spatial point r 0 is: 



£ 2 (r 0 + p) = £ 2 (r 0 )+ 



cE. 



3t 



1 T <?E 2 
P + -P — t- 



>0 



p+ 



sa+b r p + ^p T Cp 



(8) 



Where p is an incremental spatial offset to the point of interest, b is the 3-by-1 
gradient vector of E 2 at r 0 , and C is the 3-by-3 Hessian. Taking the derivative of 
this expansion with respect to p and setting it equal to zero yields: 

p = -C b (9) 
Replacing r 0 with r 0 + p, this evaluation can continue iteratively until p converges 
to (0,0,0), whence r 0 is identified as a stationary point of E 2 , and presumably a 

maximum, if the initial point was sufficiently close. 

The computation of the gradient and Hessian can be effected in terms of 
the precomputed H matrices using equations (6) and (7). The gradient at 
r = {x,y, z) is, by definition 



dE 



2 _ 



dt 



dE 2 
cE 2 
cE 2 



(10) 



where the first derivative of E 2 with respect to x is represented by 



dx dx v i 



= 2s' 



H(^r^(l-H(rH)-r) s 



(11) 



[341 1 4-8007/SL023650.065] 



-13- 



9/10/03 



and similarly for derivatives with respect to y and z. The Hessian at r = (x,y,z) 
is, by definition 



2 _ 



d*E 2 d L E 2 <?E 2 



dx dxdy dx dz 
d L E 2 c?E 2 $E 2 



dxdy dy dy dz 
d i E 2 &E 2 d L E 2 



dx dz dy dz dz 2 
The second partial with respect to x and y is 



d 2 E, 



dx dy dx dy 



s 7 H(H r H) H r s 



- * |(.-e(h'b )- r )§ (v ft )-' -h ( r h)- 1 te ) 



(12) 



(13) 



-t(ra)-^:a(rH)^i-s(r*)-r) 
-H(rH)-'^H(rH)-|l(i-H(rH)-'H^) 

- H ( H^h) ' ^(i- H ( H'H)" H r )^(H'H )" H r 

tS (4'S)'fH(H'Sr6l). 

The other entries in the Hessian can be represented similarly. 
[0035] It has been found that the computation of the gradient and Hessian can be 

effected very efficiently in terms of the precomputed H matrices, using equation 
(6). 

Weighting of the Optimization Metric 
[0036] As is known in the art of numerical optimization, when the measurements 

from some channels are more reliable than others, which may be due to 
interference, correlated noise, or nonuniform noise, then equation (2) may be 
modified according to 



£, = (Bm - s) r CT 1 (Bm - s) 



(14) 
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where C is the spatial covariance matrix of the additive noise or interference. A 
whitening matrix W can be defined such that W r W = C~\ in which case equation 
(14) can be expressed 



It is thus seen that simply replacing H with WH and s with Ws in the foregoing, 
and, in particular, in the objective function of equation (4), suffices to account for 
the weighting in the error metric in equation (14) (as described in S.M. Kay, 
Fundamentals of Statistical Signal Processing: Estimation Theory (Prentice Hall 
PTR 1993), which is hereby incorporated by reference). The whitening matrix is 
not dependent on the location of the marker, so all reference measurements are 
weighted equally. The precomputed H matrices associated with a grid of points 
could be weighted during precomputation to eliminate the need to do so at the 
time of localization. It has been found in practice, however, that the weighting can 
be effected in real time, and this is a preferred embodiment when weighting is 
needed or desired for better performance in the presence of noise or interference. 
It will be appreciated that the interpolation and maximization techniques disclosed 
above will work unchanged with weighted H matrices, and in practice the 
computation to weight the H matrices can be done either before or after 
interpolation. 

Reducing the Effects of Interference Using Column Augmentation 
[0037] In practice there are sources of magnetic interference, such as CRTs and 

motors, that can reduce the accuracy of the location system. When there is 
knowledge of the relative amplitudes of interfering signals across the sense 
channels, but not their absolute amplitudes, then the signal model of equation (1) 
may be modified according to 



where the columns of the augmentation matrix A are specified to match the 
interfering signals to within a scale factor, and the vector a contains the unknown 
scale factors associated with each column. Ideally, there is a one-to-one 
correspondence between the interfering sources and the columns of the 



E x = (WHm - Ws) r (WHm - Ws) 



(15) 



h = [H> A] 



a 



(16) 
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augmentation matrix, but in practice the columns may be chosen to correspond to 
components of interfering sources, in which case there may be more columns than 
interferes. 

[0038] It will be appreciated that a is independent of the marker location and is of 

no interest for localization. It will be further appreciated that the optimization 
procedure described above can be applied with each set of reference 
measurements H replaced by the augmented set [HA] in equation (4). 
Alternatively, if whitening matrices are used, then WH is replaced by W[HA]. In 
either event, some of the degrees of freedom of the location system can be used 
to account for the interference, hence, substantially removing its effect on the 
estimate of the marker location. 

[0039] For example, if the interfering source is a CRT displaced from the array of 

sensors by several feet, the interfering magnetic field will exhibit, to the zeroth 
order, a constant magnitude across the array. In this case, the augmentation 
matrix A consists of a single column in which each entry is a constant. The 
numerical value of the constant is unimportant to the operation of the method. A 
first order description of the CRT interference would include gradient terms across 
the array. In this case, two more columns would be added to the augmentation 
matrix A , for a total of three. Assuming the array is planar and described by a 
Cartesian coordinate system, then the first additional column could consist of 
entries proportional to the x coordinate of the associated sensor, and the second 
additional column could consist of entries proportional to the y coordinate. Once 
again, it is the proportionality, rather than the specific numerical values, that is 
pertinent. Such an approach can be extended to a second order, requiring three 
additional columns (proportional to x 2 , xy, and y 2 ), and a higher order if required. 
The inventors have found that this approach is effective in suppressing the 
detrimental effects of interfering sources that are relatively far from the array of 
sensors, for example distances of greater than two meters from an array of 40cm x 
40cm, even when the interference is of relatively large magnitude. The inventors 
have also found that this technique is effective in suppressing the detrimental 
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effects of conductors near the localization volume, which can often be modeled as 
interference. 

[0040] When it is possible to know the relative spatial characteristics of the 

interference more precisely, then the columns of the augmentation matrix can be 
specified to match, to within a scale factor, these characteristics. This 
embodiment has the advantage of suppressing interference using fewer degrees 
of freedom of the location system. However, this embodiment may be more 
sensitive to inaccuracies in the spatial characteristics of the interference, and if this 
knowledge is imperfect, then it may be preferred to use low order terms in the 
columns of the augmentation matrix as specified above. 

[0041] It will be appreciated that the augmentation columns are independent of 

marker location, so there is little processing required to augment the H matrices in 
the objective function of equation (4) prior to its actual computation. It will likewise 
be appreciated that the interpolation and maximization techniques disclosed above 
will work unchanged with augmented H matrices, and that column augmentation 
can be applied either before or after interpolation. The inventors have found that 
augmentation by one, three, or six columns (zeroth, first, or second order) can be 
effected without compromising the ability of the location system to function in real 
time. 

[0042] There is an alternative mathematical formulation to column augmentation 

that may be preferred in some embodiments. Using the signal model of equation 
(16) in the error metric of equation (14) yields 

^(Hm + Aa-s/CT'^m + Aa-s) (17) 
When this is minimized with respect to the vector a and the result substituted 
back, the metric can be written 

E x = (PWHm - PWs/ (PWHm - PWs) ( 1 8) 

where the projection matrix P is given by 

P = I - (WA) [(WA) r (WA) J' (WA) r ( 1 9) 

and I is the identity matrix. This is the projection matrix associated with the null 
space of the columns of WA. When W is associated with uniform spatially white 

[341 1 4-8007/SL023650.065] - 1 7- 9/1 0/03 



noise, then it can be replaced in equations (18) and (19) by the identity matrix. It 
will be noted that the optimization procedure described above can be applied with 
each set of reference measurements H replaced by the set PWH and the vector 
of actual measurements s with PWs in equation (4), and in particular the 
interpolation and maximization techniques can be applied unchanged. While this 
is mathematically equivalent to column augmentation, it is effected as a weighting 
of the H matrices rather than an augmentation. This approach may be preferred 
in practice as it unifies column augmentation and weighting. It will be appreciated, 
however, that this technique is distinct from equation (14) as P is singular and 
does not possess an inverse. In this regard, the method differs from those in the 
prior art of magnetic localization, as neither P nor PW can be representative of 
the inverse of a spatial noise correlation matrix. More specifically, the method 
disclosed herein is not intended to whiten the spatial noise, as in the prior art, but 
rather to suppress completely any interference that is a linear combination of the 
columns of the augmentation matrix. 
Components of the Location System 

[0043] Figure 3 illustrates components of the location system in one embodiment. 

The location system 310 receives sets of actual measurements from the generate 
actual measurements component 320. The sensor array is positioned near the 
object 340 (e.g., human body) in which the marker 350 has been implanted. In 
some embodiments, the marker may be affixed to the object, rather than implanted 
in the object. More generally, the marker is associated with the object so that the 
location of the marker is indicative of the location of the object. The bounding box 
360 indicates the volume that contains the marker. The location system includes a 
generate H's component 311, a collection of H's 312, a find bounding sub-box 
component 313, and a find marker point within bounding sub-box component 314. 

[0044] The location system may be implemented on a computer system that 

includes a central processing unit, memory, input devices (e.g., keyboard and 
pointing devices), output devices (e.g., display devices), and storage devices (e.g., 
disk drives). The memory and storage devices are computer-readable media that 
may contain instructions that implement the location system. In addition, the data 
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structures may be stored or transmitted via a data transmission medium, such as 
the Internet, a local area network, a wide area network, or a point-to-point dial-up 
connection. 

Search for Closest Grid Point 

[0045] Figure 4 is a flow diagram illustrating the processing of the find bounding 

sub-box component in one embodiment. The component identifies the bounding 
sub-box by first locating the grid point (i.e., the closest grid point) within the 
bounding box with the maximum value of the objective function. The component 
then identifies one of the eight sub-boxes adjacent to the closest grid point that is 
expected to contain the marker point. This component is passed the H matrices 
and the corresponding spatial grid spanning the bounding box. In block 401 , the 
function invokes the find closest grid point component. That component returns 
the closest grid point within the bounding box. In block 402, the component 
invokes the find adjacent sub-box component passing the closest grid point and 
the H matrices. That function returns the lower left forward point of the bounding 
sub-box. The component then completes. 

[0046] Figure 5 is a flow diagram illustrating the processing of the find closest grid 

point component in one embodiment. The component is passed the H matrices 
and the corresponding spatial grid spanning the bounding box, as well as an initial 
grid point. In this embodiment, the component implements a hill-climbing 
algorithm starting at the initial grid point. The component uses an iterative 
technique that initially sets the current grid point to the initial grid point and 
calculates the objective function at the initial grid point. The component then 
calculates the objective function at each of the six grid points adjacent to the 
current grid point (i.e., two grid points adjacent to each of the x, y, and z 
directions). The component selects the adjacent grid point whose value of the 
objective function is greatest. If the value of the objective function at that selected 
grid point is greater than the value of the objective function at the current grid 
point, the component sets the new current grid point to the selected grid point and 
repeats the process of selecting the adjacent grid point with the greatest value of 
the objective function. If, however, the value of the objective function at the 
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selected grid point is not greater than the value of the objective function at the 
current grid point, then the solution has converged on the current grid point and 
the component returns the current grid point as the closest grid point. In block 
501, the component initializes the current point (i.e., pointO) to the initial grid point. 
In block 502, the component initializes the maximum value of the objective 
function found so far (i.e., E MAX ) to the value of the objective function at the initial 
grid point and sets the maximum adjacent grid point (i.e., pointl) to the initial grid 
point. In blocks 503-513, the component loops following the grid point adjacent to 
the current grid point whose value of the objective function is greatest. In block 
503, the component sets the grid point in the x direction (i.e., pointX) to the 
adjacent grid point in the x direction with the larger value of the objective function. 
In decision block 504, if the value of the objective function at the point in the x 
direction is larger than the maximum value of the objective function encountered 
so far, then the component continues at block 505, else the component continues 
at block 506. In block 505, the component sets the maximum adjacent grid point 
to the point in the x direction and sets the maximum value of the objective function 
to the value of the objective function at the point in the x direction. In blocks 506- 
508 and 509-51 1 , the component performs similar processing for the adjacent grid 
points in the y and z directions, respectively. In decision block 512, if the current 
grid point is the same as the maximum adjacent grid point, then none of the 
adjacent grid points have a value of the objective function that is greater than that 
of the current grid point and the current grid point represents the closest grid point. 
If so, the component returns, else the component continues at block 513. In block 
513, the component sets the current grid point to the maximum adjacent grid point 
and loops to block 503 to process the six grid points adjacent to the new current 
grid point. 

[0047] Figure 6 is a flow diagram illustrating the processing of the find adjacent 

sub-box component in one embodiment. The component is passed the closest 
grid point and the H matrices. The component identifies the sub-box that contains 
the marker point. In one embodiment, the component identifies the bounding sub- 
box as the sub-box whose objective function values in the x, y, and z directions 
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are greatest. The component then calculates and returns the lower left forward 
grid point of that bounding sub-box to identify the bounding sub-box. In blocks 
601-603, the component invokes the find neighbor maximum point passing the 
closest grid point, an indication of one of the three dimensions, and the H 
matrices. These components return an indication of the grid points in each of the 
dimensions that have the largest value for their objective functions. In block 604, 
the component calculates the lower left forward point of the bounding sub-box and 
returns. 

[0048] Figure 7 is a flow diagram illustrating the processing of the find neighbor 

maximum point component in one embodiment. The component is passed the 
closest grid point, a direction, and the H matrices. In block 701 , the component 
identifies the coordinates of the neighbor grid point in the lower direction. For 
example, if the closest grid point is (10,10,10) and the direction is the x direction, 
then the identified coordinates are (9,10,10). In block 702, the component 
evaluates the objective function for the lower neighbor grid point. In block 703, the 
component identifies the upper neighbor grid point. In block 704, the component 
evaluates the objective function for the upper neighbor grid point. In decision 
block 705, if the value of the objective function for the lower neighbor grid point is 
greater than the value of the objective function for the upper neighbor grid point, 
the component returns the lower neighbor grid point, else the component returns 
the upper neighbor grid point. 
Search for Marker Point 

[0049] Once the bounding sub-box is identified, the location system calculates the 

marker point using interpolated H's within the bounding sub-box based on the H's 
of the grid points. The location system uses a Taylor series expansion of the 
objective function based on the interpolated H's to converge on the off-grid point 
with the maximum value of the objective function. The location system uses an 
iterative technique that starts at an initial point within the bounding sub-box and 
calculates spatial offsets according to equation (9) until the spatial offsets 
converge to zero at a point which represents the marker location. 
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[0050] Figure 9 is a flow diagram illustrating the processing of the find marker point 

within a bounding sub-box component in one embodiment. The component is 
passed the H matrices, the spatial grid, the bounding sub-box, and an initial point. 
The component initializes the current point to the initial point. The component then 
calculates the spatial offset of the current point according to equation (9). If the 
spatial offset is not within a certain tolerance of a zero offset, then the component 
sets the new current point to the previous current point plus the spatial offset. 
Otherwise, the component has converged to a solution and returns the current 
point as the marker point. In block 901 , the component sets the current point (i.e., 
pointO) to the initial point. In block 902, the component calculates the gradient of 
the objective function at the current point according to equation (1 1 ). In block 903, 
the component calculates the Hessian of the objective function at the current point 
according to equation (13). In block 904, the component calculates the 
incremental spatial offset according to equation (9). In decision block 905, if the 
spatial offset is within a certain tolerance of zero, then the component returns the 
current point as the marker pointer, else the component continues at block 906. In 
block 906, the component sets the new current point to the previous current point 
plus the spatial offset. The component then loops to block 902 to calculate the 
next spatial offset. 

(0051] One skilled in the art will appreciate that although specific embodiments of 

the location system have been described herein for purposes of illustration, 
various modifications may be made without deviating from the spirit and scope of 
the invention. Accordingly, the invention is not limited except by the appended 
claims. 



[341 1 4-8007/SL023650.065] 



-22- 



9/10/03 



